Strain versus stress in a model granular material: a Devil's staircase. 
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The series of equilibrium states reached by disordered packings of rigid, frictionless discs in two 
dimensions, under gradually varying stress, are studied by numerical simulations. Statistical proper- 
ties of trajectories in configuration space are found to be independent of specific assumptions ruling 
granular dynamics, and determined by geometry only. A monotonic increase in some macroscopic 
loading parameter causes a discrete sequence of rearrangements. For a biaxial compression, we show 
that, due to the statistical importance of such events of large magnitudes, the dependence of the 
resulting strain on stress direction is a Levy flight in the thermodynamic limit. 

PACS numbers: 83.70.Fn,05.40.-a,45.05.-|-x 



The mechanical properties of granular media are cur- 
rently an active field of research, both in the condensed 
matter physics and in the mechanics and engineering 
communities P, 0, 01- 

Granular packings close to mechanical equilibrium are 
traditionnally modelled, in the framework of continuum 
mechanics, with elastoplastic constitutive laws 0, 
i.e., incremental stress-strain relations. Such laws, de- 
spite their practical success, were never clearly related 
to grain-level mechanics. Moreover, cohesionless gran- 
ular systems seem to be quite different from ordinary 
solids. Many experimental, theoretical and numerical 
studies d, 0, d, 0, Hi have recently been devoted to the 
peculiar features of stress transmission in granular sys- 
tems at equilibrium, with correlations over length scales 
significantly larger than the grain size. 

Observations of displacement fields and strains, as 
the system moves from one equilibrium to another, are 
scarcer. Systems of rigid grains are expected to deform 
because of rearrangements of the packing, rather than 
contact elasticity. How such rearrangements average to 
produce a macroscopic strain, related to stress variations, 
remains rather mysterious. The rather singular, unilat- 
eral form of the local interaction in such systems led some 
authors Q to expect quite unusual macroscopic proper- 
ties, for which the very concept of strain, so familiar in 
mechanics of solids, would be irrelevant. 

Direct grain-level approaches are, in principle, possible 
by numerical simulations. However, one has then to de- 
fine a complete mechanical model to enable a calculation 
of particle trajectories. In practice, dynamical parame- 
ters ruling energy dissipation are often chosen according 
to computational convenience as much as physical accu- 
racy. It would be desirable to assess the influence of such 
choices on the results. 

The present numerical study addresses those problems, 
as follows. 

Disordered, dense assemblies of rigid, circular, friction- 
less discs, are prepared by isotropic compaction. The 
force law reduces to the condition that contacts transmit 
repulsive normal forces of unknown magnitude. Then, 
the direction of the load is gradually altered, thus simu- 



lating the biaxial test of fig. [T] Exploiting the isostaticity 
property [13, El, of such systems, which is exactly 
satisfied provided impenetrability is enforced accurately 
enough, we designed a prescription for the compu- 
tation of sequences of equilibrium states, that we call 
the geometric quasi-static method (GQSM), in which the 
only inputs are the geometric data. The strain versus 
stress evolution in biaxial compression tests is recorded, 
and a statistical analysis of fluctuations and system size 
dependence is carried out. Then the GQSM predictions 
are compared with those of a standard molecular dynam- 
ics (MD) method. 

First, samples of different sizes are generated (51 sam- 
ples of N=1024 discs, 23 with N=1936, 10 with N = 3025 
and 15 with N=4900). Disc diameters are uniformly dis- 
tributed between 0.5 and 1 (the largest diameter is cho- 
sen as unit length), and packed in a loose state within 
a rigid, square box. After some amount of random mix- 
ing (using some dynamical method with energy conserva- 
tion) we proceed to the isotropic compaction: two of the 
walls, 1 and 2 on fig. [1] are now mobile, and submitted to 
compressive external forces -Fi (along the x axis on the 
figure) and F2 — Fi (along the y axis). Stress compo- 
nents (Til = F1/L2 and (722 = F2/L1 are kept constant, 
equal to p, while the system lengths along directions 1 
and 2 {Li and L2) decrease. To produce a dense, iso- 
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FIG. 1: Sketch of a biaxial compression test. 



static equilibrium state we use the 'lubricated granular 
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dynamics' method of refs. 0, [l3|- Then, the initial, ref- 
erence configuration of the biaxial experiment is ready: 
F2 ~ {p + q)Li is gradually increased while Fi ~ 
stays constant, strain components are defined as the rela- 
tive decrease of lengths (Li)i<i<2, with reference to their 
initial values (L°)i<i<2 as en = — ALi/L°. One also de- 
fines the volumetric jisj strain as = —en — 622- We 
use units such that p = 1. 

Our essential result here is the obtention of the e{q) 
curves in the thermodynamic limit, as loading param- 
eter q increases monotonically, at constant p. Let us 
first describe the GQSM procedure. One starts from an 
equilibrium state in which the force-carrying structure 
is isostatic. This means that the equilibrium conditions 
are sufficient to compute all contact force values, on the 
one hand (the structure is devoid of hyperstaticity, it is 
not 'overbraced' or 'overconstrained' [lO|), and that the 
force-carrying structure is rigid (devoid of mechanisms or 
'floppy modes' [13]) on the other hand. The first of these 
two properties stems from the condition that two grains 
need be exactly in contact to transmit a force to one an- 
other [111, and cannot interpenetrate. It entails that 
force values, once equilibrium positions are known, are 
geometrically determined, all material properties being 
irrelevant in the limit of rigid grains. The second prop- 
erty is satisfied, for stability reasons, because the grains 
are circular and contacts do not withstand tension [l^ . 
It entails that an assembly of rigid discs will not deform 
at all until some initially active contact opens. This can- 
not occur as long as contact forces are compressive, since 
this would require the potential energy to increase from 
equilibrium. As soon as one contact force vanishes, this 
contact will open [13], because the resulting motion cor- 
responds to an instability. Hence the following algorithm: 

(1) At equilibrium, as q increases from its initial value 
qo, the contact forces depend linearly on q (equilibrium 
equations are linear). When q reaches some value qo + Sq, 
the force vanishes in one contact, say Iq. 

(2) Open Iq, all other contacts being maintained. Due 
to isostaticity, this entirely determines the initial direc- 
tion of motion for the whole structure. Keep moving the 
grains with the same prescription. 

(3) When another contact, say ^i, closes, the new con- 
tact structure (the old one, minus Iq, plus h) is isostatic 
and may carry the load with geometrically determined 
contact forces. If there is no traction, a new equilibrium 
state has been found: go back to step 1. Otherwise, pick 
up the largest traction, call the corresponding contact Iq 
and go back to step 2. 

This procedure determines a series of equilibrium con- 
figurations, that are separated by rearrangements occur- 
ring for discrete values of q. The strain versus stress curve 
is a staircase (see fig. [2]). As long as the same contacts 
carry the load, the system does not deform; as soon as 
a rearranging event occurs, strain variables jump to the 
values corresponding to the next equilibrium configura- 
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FIG. 2: 'Axial' strain £22 versus deviator q in one sample with 
N=4900. 



tion. This algorithm clearly involves, in steps 2 and 3, an 
arbitrary ingredient: the prescription that contacts open 
one by one. The main merit of GQSM, however, is that 
it does not introduce parameters other than geometric 
ones. 

We now focus on the rise of 622 with q, close to the ori- 
gin, and ask whether the staircase approaches a smooth 
curve in the thermodynamic limit. (Its initial slope, if 
finite, would be the effective compliance of the material). 
To do so, one studies the statistics of stress (Sq) and 
strain ((5e22, ^eii) steps. 

Successive Sq and Se values are found to be indepen- 
dent, and the width Sq of a stability interval is not cor- 
related to the following strain steps. Throughout the 
investigated q interval, the probability distributions of 
increments 5q, 6622, and 6W = pSey — qSe22, which is the 
variation in potential energy corresponding to the cur- 
rent load, do not appreciably change [l^ . No significant 
difference between samples is observed either. 

Both q and 622 values reached at a given stage can thus 
be regarded as sums of equidistributed independent ran- 
dom increments. The distribution of stress increments 
dq is displayed on fig. [31 It decays exponentially for 
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FIG. 3; Probability distribution function for rescaled stress 
increments 5g(A'^/1024)^'^'^ for the 4 values of system size A*'. 
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large Sq, and is shifted to smaller and smaller values as 
N increases, so that the probability distribution of the 
rescaled increment SqN" is size-independent. We denote 
as SqQ its average. The exponent can be estimated as 
a — 1.16 ± 0.04. Stability intervals shrink to zero as N 
increases: in the thermodynamic limit, any macroscopic 
load variation entails some motion of the grains. This 
property of packings of frictionless rigid grains, known as 
fragility [1, [l3| , is unambiguously established by our sim- 
ulations. The value of a should be related to the shape 
of the force distribution for small values, and to the vary- 
ing sensitivity of the contact forces to macroscopic load 
increments. The classical central-limit theorem, applied 
to q increments, will relate for large systems the number 
of steps M from the beginning of the biaxial compression 
to the current value of q, as 

M - -^iV". (1) 

Oqo 

The distribution of strain increments Se22 is shown on 
figure m A reduced variable 6622/ Sen can also be de- 




FIG. 4: Probability distribution function for rescaled strain 
increments 5e22(A''/1024)^'^® for the 4 values of TV. The insert 
displays the probability distribution function of the number c 
of lost contacts in one strain step for each system size. 

fined, with an iV-independent probability distribution. 
Remarkably, the (power-law distributed) number c of 
contact losses (or gains) in one strain step, which tends to 
grow with Se22, does not show any significant size depen- 
dence, and contacts that simultaneously open or close are 
homogeneously scattered throughout the sample, what- 
ever N. We estimated /3as/3 = 2.18±0.12. Unhke for 
stress increments, the distribution now decays as a power 
law: 

p{Se22) - {Se22)-^^+''\ 

with ^ = 0.46 ± 0.03. As /i < 1, it should be remarked 
that, although the typical strain increment decreases as 
N^^ in the limit of large samples, the average strain in- 
crement does not exist. The standard central limit the- 
orem does not apply, but, due to the power law decay of 



the distribution function, one may resort to a generalized 
central limit theorem [16{ . In our case, for large numbers 
M of increments, the asymptotic form of their sum £22 is 

£22 5toN-^M^^, (2) 

with an Af-independent random variable ^ abiding by an 
asymmetric Levy distribution of index ^. 

The stress-strain relationship is obtained on combining 
eqns [T] and [D 

£22 = ^q^N-^+H- (3) 
(5go)" 

The presence of ^ in relation[3]implies that the strain ver- 
sus stress curve will never express a deterministic depen- 
dence. Some size effect - a dependence on iV - remains 
in the thermodynamic limit, unless one has: 

a = P^Ji (4) 

Our estimates of a, (3 and /i are compatible with this rela- 
tion. If it is satisfied, then the distribution of axial strain 
increments Ae22 corresponding to a given, fixed, q incre- 
ment, Ag, should no longer depend on the system size 
as soon as N is large enough for Ag to involve, typically, 
sufficiently many elementary rearranging steps. This was 
checked, confirming (see fig. [5|), within statistical errors, 
relation m and the absence of size effects. 

In order to compare the GQSM results to the pre- 
dictions of more conventional methods, we also simu- 
lated biaxial compressions on samples of 1024 and 3025 
discs using MD, successively imposing deviator incre- 
ments Ag — 10~^. Contacts obey an elastic unilateral 
law, with a normal stiffness equal to 10^. After each 
stress step, one waits for a new equilibrium state, request- 
ing forces on each grain to balance up to an accuracy of 
10~^. Successive strain increments Aef^^ are uncorre- 
lated and distributed according to the same probability 
law, which coincides (within statistical uncertainties) as 
shown on fig. O with that of increments Ae22 obtained 
for the same value of Ag by GQSM. Similarly, the(power- 
law) distribution of the fraction of the total number of 
contacts that change within one fixed increment Ag does 
not depend on TV (fig. [5]). We therefore conclude that, 
in a biaxial test, the axial strain £22 dependence on de- 
viator g, as expressed by relations [3] and HI is a Levy 
stochastic process. It does not become deterministic in 
the thermodynamic limit, but it is devoid of size effect: 
the strain versus stress curve approaches a Devil's stair- 
case with a dense set, on the g axis, of discontinuities of 
random magnitudes. Moreover, it appears not to depend 
on dynamical parameters (at least within the accuracy 
of the present study) and to be determined by the sole 
system geometry. These results apply, without appre- 
ciable change in the statistics, throughout the interval 
< g/p < 0.2. 
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FIG. 5: Probability distribution function for strain increments 
Ae22 corresponding to Ag = 10~^p obtained for the 4 values 
of by GQSM (continuous curves) and for N=1024 (crosses) 
and N=3025 (square dots) by MD. The insert shows the dis- 
tribution of C'/N, where C is the number of lost contacts in 
one increment Ag, for the different system sizes. 



The evolution of (which can be of either sign, typ- 
ically one order of magnitude smaller than £22) is some- 
what more complicated, since successive increments Sey 
are not equidistributed (unlike SW). Leaving a discus- 
sion of volumetric strains to a future publication, let us 
further comment here the behavior of axial ^is'l strains. 

Our results preclude the existence of a constitutive law 
in the usual sense. As q is increased, £22 is typically of or- 
der g^/'', but its actual value is essentially impredictible, 
with the remarkable consequence that macroscopic mod- 
els for granular mechanics should be of a stochastic, 
rather than deterministic, nature. Of course, such a be- 
havior might be limited to frictionless grains, and is thus 
perhaps more relevant for some colloidal glasses than for 
sand, although important fluctuations and high noise lev- 
els were sometimes reported for soil materials or glass 
beads. One does obtain an extremely noisy curve on 
numerically compressing our system at constant strain 
rate instead of controlling the stress: understandably, 
between two equilibrium states that can be relatively far 
apart, the contact structure docs not consistently oppose 
any given level of stress, and the observed response de- 
pends on dynamical parameters. Laboratory tests, at 
constant strain rate, on such systems, would be sensi- 
tively influenced by the apparatus itself. It will be inter- 
esting to study the dependence of parameter fi on grain 
shape and polydispersity. Preliminary MD results on 
dense random packings of monodisperse spheres in 3D 
yielded Levy-distributed large strain steps with a value 
of ^ in the 0.4-0.6 range. 

Physically, large strain increments are due to re- 
arrangements involving many contact changes (GQSM 
steps 2 and 3 have to be repeated). Had we resorted 



to the approximation of small displacements (ASD), in 
which all relevant quantities are dealt with to leading 
order in the displacements, then, as shown in ref. 
no iteration of steps 2 and 3 would have been necessary. 
Within the ASD, contacts are replaced one by one, and 
we have checked that the resulting strain increment dis- 
tribution does possess an average value. This geometric 
approximation would ensure uniqueness of the equi- 
librium state and a one-to-one correspondence between 
stress and strain. The ASD was used for slightly poly- 
disperse discs on a triangular lattice, in which case it can 



be justified and constitutive laws are obtained [ij, |14| 
(the second ref. [lT| proposes an independent implemen- 
tation of the GQSM, with the ASD). The distinctive me- 
chanical features of model granular assemblies that are 
reported here are thus due to the motion of the system in 
configuration space along a tortuous path between local 
minima in a complex potential energy lanscape. Coop- 
erative rearranging events are also observed in glassy re- 
laxation 17[, which might open interesting perspectives. 

We are grateful to Jean-Philippe Bouchaud for useful 
suggestions. 
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